Competition-induced greater trait variability obscures the trait-growth relationship

1 Core function

# c('#304D30', '#4F6F52', '#739072', '#86A789', '#D2E3C8', '#B6C4B6',
# '#DFD0B8', '#DEAA79', '#AF8F6F', '#6F4E37')
calc_annual_trait_RGR <- function(data, trait, RGR = "RGR") {
    result <- data.frame()
    for (compete in unique(data$Compete)) {
        for (year in unique(data$comp_year)) {
            subset_data <- data %>%
                filter(Compete == !!compete, comp_year == !!year) %>%
                select(Nenv, sp, ring.num, Compete, growth_time, Light, Moisture,
                  AP, trait, RGR, paste0(trait, "_minmax_scaled")) %>%
                na.omit()

            if (nrow(subset_data) > 5) {
                # Fit the mixed-effects model
                model_formula <- as.formula(paste0(RGR, " ~ ", trait, " *(Light + Moisture + AP) + (1 | sp) + (1|Nenv)"))
                model <- lmer(model_formula, data = subset_data)
                # Extract slope (first fixed effect coefficient)
                slope <- fixef(model)[2]
                TLslope <- fixef(model)[6]
                TMslope <- fixef(model)[7]
                TPslope <- fixef(model)[8]

                # Calculate marginal R² using `MuMIn` package for lmer models
                library(MuMIn)
                r2m <- r.squaredGLMM(model)[, "R2m"]
                r2c <- r.squaredGLMM(model)[, "R2c"]

                sp_slopes <- c()  # Initialize a vector to store slopes for each species
                # Iterate through each species
                for (species in unique(subset_data$sp)) {
                  # Subset data for the current species
                  species_data <- subset_data %>%
                    filter(sp == species) %>%
                    na.omit()
                  # Ensure there is enough data to fit the model
                  if (nrow(species_data) > 5) {
                    # Define the model formula dynamically for the given trait
                    # and RGR
                    model_formula <- as.formula(paste0(RGR, " ~ ", trait))
                    # Fit the linear mixed-effects model for the current
                    # species
                    sp_model <- summary(lm(model_formula, data = species_data))
                    # Extract the slope for the trait (first fixed-effect
                    # coefficient)
                    sp_slope <- sp_model$coefficients[2, 1]
                    # Append the slope to the slopes vector
                    sp_slopes <- c(sp_slopes, sp_slope)
                  }
                }

                # Calculate the mean and standard deviation of slopes
                normalized_slopes <- (sp_slopes - min(sp_slopes))/(max(sp_slopes) -
                  min(sp_slopes))
                mean_slope <- mean(normalized_slopes, na.rm = TRUE)
                sd_slope <- sd(normalized_slopes, na.rm = TRUE)
                # Calculate the coefficient of variation (CV) of slopes
                cv_slope <- sd_slope/mean_slope
                # Extract p-value for the main effect of the trait
                p_value <- summary(model)$coefficients[2, "Pr(>|t|)"]

                com_itv <- CV4(subset_data[, paste0(trait, "_minmax_scaled")])
                # sp_itv <- mean(tapply(subset_data[, paste0(trait,
                # '_minmax_scaled')], list(subset_data$sp), CV4), na.rm = T)

                subset_data$new_com <- "Com"
                kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$new_com)  # Kernel density estimates of trait for individuals
                functional.indices.trait <- REND(TPDs = kernel.trait)

                sp.kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$sp)  # Kernel density estimates of trait for individuals
                sp.fun.indices <- REND(TPDs = sp.kernel.trait)
                sp_itv <- mean(sp.fun.indices$species$FRichness)
                ol <- dissim(sp.kernel.trait)
                dis <- ol$populations$dissimilarity
                inter_dis <- mean(dis[upper.tri(dis) | lower.tri(dis)])

                mean_trait <- mean(subset_data[, paste0(trait, "_minmax_scaled")],
                  na.rm = T)
                sp_meantrait <- tapply(subset_data[, paste0(trait, "_minmax_scaled")],
                  list(subset_data$sp), mean, na.rm = T)
                # inter_dis <- CV4(sp_meantrait)

                RGR_var <- sd(subset_data$RGR, na.rm = T)
                RGR_mean <- mean(subset_data$RGR, na.rm = T)

                result <- rbind(result, data.frame(Compete = compete, sample_size = nrow(subset_data),
                  trait = trait, comp_year = year, estimate = slope, TLestimate = TLslope,
                  TMestimate = TMslope, TPestimate = TPslope, TYestimate = NA, Marginal_R2 = r2m,
                  Conditional_R2 = r2c, p_value = p_value, com_ITV = com_itv, sp_ITV = sp_itv,
                  Fric = functional.indices.trait$species$FRichness, Fdiv = functional.indices.trait$species$FDivergence,
                  mean_trait = mean_trait, inter_dis = inter_dis, RGR_var = RGR_var,
                  RGR_mean = RGR_mean, sp_slope_cv = cv_slope))
            }
        }
    }
    return(result)
}

calc_trait_RGR <- function(data, trait, RGR = "RGR") {
    result <- data.frame()
    for (compete in unique(data$Compete)) {
        subset_data <- data %>%
            filter(Compete == !!compete) %>%
            select(Nenv, sp, Compete, growth_time, Light, Moisture, AP, trait, RGR,
                paste0(trait, "_minmax_scaled")) %>%
            na.omit()

        subset_data$growth_time <- scale(as.numeric(subset_data$growth_time))

        if (nrow(subset_data) > 2) {
            # Fit the mixed-effects model
            model_formula <- as.formula(paste0(RGR, " ~ ", trait, " *(Light + Moisture + AP + growth_time) + (1 | sp) +(1|Nenv)"))
            model <- lmer(model_formula, data = subset_data)
            # Extract slope (first fixed effect coefficient)
            slope <- fixef(model)[2]
            TLslope <- fixef(model)[7]
            TMslope <- fixef(model)[8]
            TPslope <- fixef(model)[9]
            TYslope <- fixef(model)[9]

            # Calculate marginal R² using `MuMIn` package for lmer models
            library(MuMIn)
            r2m <- r.squaredGLMM(model)[, "R2m"]
            r2c <- r.squaredGLMM(model)[, "R2c"]

            sp_slopes <- c()  # Initialize a vector to store slopes for each species
            # Iterate through each species
            for (species in unique(subset_data$sp)) {
                # Subset data for the current species
                species_data <- subset_data %>%
                  filter(sp == species) %>%
                  na.omit()
                # Ensure there is enough data to fit the model
                if (nrow(species_data) > 5) {
                  # Define the model formula dynamically for the given trait
                  # and RGR
                  model_formula <- as.formula(paste0(RGR, " ~ ", trait))
                  # Fit the linear mixed-effects model for the current species
                  sp_model <- summary(lm(model_formula, data = species_data))
                  # Extract the slope for the trait (first fixed-effect
                  # coefficient)
                  sp_slope <- sp_model$coefficients[2, 1]
                  # Append the slope to the slopes vector
                  sp_slopes <- c(sp_slopes, sp_slope)
                }
            }

            # Calculate the mean and standard deviation of slopes
            normalized_slopes <- (sp_slopes - min(sp_slopes))/(max(sp_slopes) - min(sp_slopes))
            mean_slope <- mean(normalized_slopes, na.rm = TRUE)
            sd_slope <- sd(normalized_slopes, na.rm = TRUE)
            # Calculate the coefficient of variation (CV) of slopes
            cv_slope <- sd_slope/mean_slope
            # Extract p-value for the main effect of the trait
            p_value <- summary(model)$coefficients[2, "Pr(>|t|)"]

            com_itv <- CV4(subset_data[, paste0(trait, "_minmax_scaled")])
            # sp_itv <- mean(tapply(subset_data[, paste0(trait,
            # '_minmax_scaled')], list(subset_data$sp), CV4), na.rm = T)
            subset_data$new_com <- "Com"
            kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$new_com)  # Kernel density estimates of trait for individuals
            functional.indices.trait <- REND(TPDs = kernel.trait)

            sp.kernel.trait <- TPDs(traits = subset_data[, trait], species = subset_data$sp)
            sp.fun.indices <- REND(TPDs = sp.kernel.trait)
            sp_itv <- mean(sp.fun.indices$species$FRichness)
            ol <- dissim(sp.kernel.trait)
            dis <- ol$populations$dissimilarity
            inter_dis <- mean(dis[upper.tri(dis) | lower.tri(dis)])

            mean_trait <- mean(subset_data[, paste0(trait, "_minmax_scaled")], na.rm = T)
            sp_meantrait <- tapply(subset_data[, paste0(trait, "_minmax_scaled")],
                list(subset_data$sp), mean, na.rm = T)
            # inter_dis <- CV4(sp_meantrait)
            RGR_var <- sd(subset_data$RGR, na.rm = T)
            RGR_mean <- mean(subset_data$RGR, na.rm = T)

            result <- rbind(result, data.frame(Compete = compete, sample_size = nrow(subset_data),
                trait = trait, comp_year = "All", estimate = slope, TLestimate = TLslope,
                TMestimate = TMslope, TPestimate = TPslope, TYestimate = TYslope,
                Marginal_R2 = r2m, Conditional_R2 = r2c, p_value = p_value, com_ITV = com_itv,
                sp_ITV = sp_itv, Fric = functional.indices.trait$species$FRichness,
                Fdiv = functional.indices.trait$species$FDivergence, mean_trait = mean_trait,
                inter_dis = inter_dis, RGR_var = RGR_var, RGR_mean = RGR_mean, sp_slope_cv = cv_slope))
        }
    }
    return(result)
}

# Core function
CV4 <- function(trait_sample) {
    trait_sample <- trait_sample[!is.na(trait_sample)]
    if (length(trait_sample) > 1) {
        N <- length(trait_sample)
        # calcualte CV^2
        y_bar <- mean(trait_sample)
        s2_hat <- var(trait_sample)
        cv_2 <- s2_hat/y_bar^2
        cv_1 <- sqrt(cv_2)
        gamma_1 <- sum(((trait_sample - y_bar)/s2_hat^0.5)^3)/N
        gamma_2 <- sum(((trait_sample - y_bar)/s2_hat^0.5)^4)/N
        bias2 <- cv_1^3/N - cv_1/4/N - cv_1^2 * gamma_1/2/N - cv_1 * gamma_2/8/N
        cv4 <- cv_1 - bias2
    } else {
        cv4 <- NA
    }
    return(cv4)
}


Comp_PCA <- function(alltraits_standard_tda, pr, trait.names, p.title, x.lim, y.lim,
    point_shape) {
    # Point data for PCA plot
    point.pca <- alltraits_standard_tda[, c("Nenv", "sp", "growth_time", "Compete",
        "PC1", "PC2")]

    # Axis data for PCA plot
    x <- "V1"
    y <- "V2"
    par.all.loading <- data.frame(Trait = trait.names, as.data.frame(matrix(pr$loadings,
        nrow = length(trait.names))))
    axis.pca <- transform(par.all.loading, v1 = 5 * (get(x)), v2 = 5 * (get(y)))
    axis.pca <- data.frame(axis.pca)

    # Calculate explained variance for PC1 and PC2
    pc1_explained <- round(100 * pr$sd[1]^2/sum(pr$sd^2), 1)
    pc2_explained <- round(100 * pr$sd[2]^2/sum(pr$sd^2), 1)

    # Create the PCA plot
    point.pca1 <- point.pca
    point.pca1$cg <- paste0(point.pca1$Compete, point.pca1$growth_time)
    point.pca1$cg <- factor(point.pca1$cg, levels = c("alone360", "alone720", "alone1080",
        "inter360", "inter720", "inter1080"))
    point.pca1$spg <- paste0(point.pca1$sp, point.pca1$growth_time)
    point.pca1$spg <- factor(point.pca1$spg, levels = c(unique(point.pca1$spg)))
    # point.pca1$compete <- factor(point.pca1$Compete, levels = c('alone',
    # 'inter'))

    ## 绘图
    p <- ggplot(point.pca1, aes(x = PC1, y = PC2)) + scale_shape_manual(values = point_shape) +
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ geom_point(aes(color
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ sp,
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ shape
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ Compete),
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ size
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 1.5,
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ alpha
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 0.7)
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ +
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ #
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ stat_ellipse(aes(group=sp,color
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ sp),fill=NA,geom='polygon'
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ ,type='norm',alpha=0.2,level
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ 0.95,show.legend
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ =
        geom_point(aes(color = sp, shape = Compete), size = 1.5, alpha = 0.7) + # stat_ellipse(aes(group=sp,color = sp),fill=NA,geom='polygon' ,type='norm',alpha=0.2,level = 0.95,show.legend = T,size=1)+ T,size=1)+
    ggtitle(p.title) + xlab(paste0("PC1 (", pc1_explained, "%)")) + ylab(paste0("PC2 (",
        pc2_explained, "%)")) + geom_segment(data = axis.pca, aes(x = 0, y = 0, xend = v1,
        yend = v2), arrow = arrow(length = unit(0.2, "cm")), alpha = 0.75, color = "black",
        linetype = 2) + geom_text(data = axis.pca, aes(v1 + 0.1, v2 + 0.1, label = Trait),
        size = 4, nudge_x = +0.2, nudge_y = +0.3) + theme(legend.text = element_text(face = "italic"),
        legend.position = "none") + theme_Publication() + scale_colour_Publication() +
        scale_fill_Publication() + xlim(x.lim) + ylim(y.lim)

    # Return the PCA plot
    return(p)
}

2 Data pre-processing

The initial dataset comprised 5188 seedling records, including 1077 dead seedlings. To ensure data quality, individuals with unrealistic growth values were removed. Specifically, seedlings with negative growth due to stem breakage or other causes, as well as those with growth rates exceeding four times the standard deviation, were excluded, resulting in the removal of 41 individuals. After this initial filtering, 4070 seedlings with valid relative growth rate (RGR) data remained.

However, during trait measurement, some seedlings lacked leaf or fine root trait data. To maintain consistency across datasets, these seedlings (a total of 740 individuals) were further excluded. Consequently, the final dataset consisted of 3330 individuals with complete functional trait and growth data. This dataset included 1952 seedlings harvested in the second year and 1378 seedlings harvested in the third year.

# Load processed growth data
load("01_data/processed_data/trait_growth_data.Rdata")

# Filter data for the 2nd and 3rd comp_years and rename RGR column
tda <- tda[tda$comp_year %in% c("2nd", "3rd"), ] %>%
  rename(RGR = RGR_H_mean) # 5188 # uses the annual average height growth rate (RGR) as an indicator of plant growth

# Count individuals with memo "S" and remove them
nrow(tda[tda$memo == "S", ]) # 1077

## [1] 1077

tda <- subset(tda, memo != "S") # 4111

# Replace negative growth values with NA and count replacements
sum(sapply(tda["RGR"], function(col) sum(col <= 0, na.rm = TRUE))) # 26

## [1] 26

tda["RGR"] <- lapply(tda["RGR"], function(col) ifelse(col <= 0, NA, col))

# Calculate mean and standard deviation for RGR
mean_rgr <- mean(tda$RGR, na.rm = TRUE)
sd_rgr <- sd(tda$RGR, na.rm = TRUE)
# Identify and count outliers beyond mean ± 4 * SD
sum(tda$RGR > (mean_rgr + 4 * sd_rgr) | tda$RGR < (mean_rgr - 4 * sd_rgr), na.rm = TRUE)

## [1] 3

# Set outlier values to NA
tda$RGR <- ifelse(tda$RGR > (mean_rgr + 4 * sd_rgr) | tda$RGR < (mean_rgr - 4 * sd_rgr), NA, tda$RGR)
# Remove individuals with NA RGR or RGR <= 0
tda <- tda[!is.na(tda$RGR) & tda$RGR > 0, ]


# Extract individuals with all valid traits (non-missing and non-zero values)
tda[trait.names] <- lapply(tda[trait.names], as.numeric) # Convert trait names to numeric
indtraits_tda <- tda %>%
  # Select columns
  select(sp, sp_abbr, Nenv, ring.num, Light, Moisture, AP, Compete, growth_time, comp_year, BM, RGR, all_of(trait.names)) %>%
  mutate(across(
    all_of(trait.names[!trait.names %in% c("LTO", "LDMC")]),
    ~ scale(log10(.))
  )) %>%
  mutate(across(
    all_of(c("LTO", "LDMC")), # "RD" is more normal
    ~ scale((.))
  )) %>%
  mutate(across(all_of(trait.names),
    ~ (. - min(., na.rm = TRUE)) / (max(., na.rm = TRUE) - min(., na.rm = TRUE)),
    .names = "{.col}_minmax_scaled"
  )) %>%
  na.omit() ## 3330

indtraits_tda$RGR <- scale(indtraits_tda$RGR)


# Replace outliers in SLA and LMA with NA
indtraits_tda$SLA[indtraits_tda$SLA < -4] <- NA
indtraits_tda$LMA[indtraits_tda$LMA > 4] <- NA
indtraits_tda$RD[indtraits_tda$RD < -5] <- NA

# Calculate effective sample sizes (non-NA and non-zero values) grouped by year and competition treatment
effective_sample_sizes <- indtraits_tda %>%
  group_by(comp_year, Compete) %>%
  summarise(across(
    all_of(trait.names),
    ~ sum(!is.na(.) & . != 0),
    .names = "{.col}"
  )) %>%
  pivot_longer(
    cols = -c(comp_year, Compete),
    names_to = "Trait",
    values_to = "Sample_Size"
  ) %>%
  unite("Time_Compete", comp_year, Compete, sep = "_") %>%
  pivot_wider(
    names_from = Time_Compete,
    values_from = Sample_Size
  )
kable(effective_sample_sizes)

Trait	2nd_alone	2nd_inter	3rd_alone	3rd_inter
LA	1008	922	692	684
SLA	1008	922	692	681
LMA	1008	922	692	681
LDMC	1008	922	692	684
LTO	1008	922	692	684
SSD	1008	922	692	684
SRL	1008	922	692	684
RD	1008	922	687	680
RTD	1008	922	692	684
SRA	1008	922	692	684

write.csv(effective_sample_sizes, file = "02_output_results/effective_sample_sizes.csv")

3 Single-trait growth models

all_indTrait_RGR_results <- data.frame()
for (trait in trait.names) {
    trait_result <- calc_trait_RGR(indtraits_tda, trait, "RGR")
    all_indTrait_RGR_results <- rbind(all_indTrait_RGR_results, trait_result)
}

indTrait_RGR_results <- data.frame()
for (trait in trait.names) {
    trait_annual_result <- calc_annual_trait_RGR(indtraits_tda, trait, "RGR")
    indTrait_RGR_results <- rbind(indTrait_RGR_results, trait_annual_result)
}

indTrait_RGR_results$comp_year <- factor(indTrait_RGR_results$comp_year, levels = c("2nd",
    "3rd"))
indTrait_RGR_results$trait <- factor(indTrait_RGR_results$trait, levels = trait.names)

all_indTrait_RGR_results$comp_year <- factor(all_indTrait_RGR_results$comp_year,
    levels = c("All"))
all_indTrait_RGR_results$trait <- factor(all_indTrait_RGR_results$trait, levels = trait.names)

4 Figure 2 R2 and slope of growth models

### Compare competition-free and competition
### Marginal R2
create_ind_trait_plot <- function(indTrait_RGR_results, y, ylab, xlab, title) {
  library(ggplot2)
  library(ggpubr) # For stat_compare_means
  library(patchwork) # If you are using patchwork for layout
  library(ggthemes) # For theme_Publication

  ggplot(indTrait_RGR_results, aes(x = Compete, y = !!sym(y), group = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.6), aes(group = Compete, linetype = Compete), width = 0.6, alpha = 1, size = 0.7) +
    geom_jitter(position = position_dodge(width = 0.6), aes(col = trait, group = Compete, shape = Compete), alpha = 0.3, size = 2.5) +
    geom_line(aes(group = trait, color = trait), size = 0.8, alpha = 0.3) +
    scale_fill_manual(values = c("#cc340c", "#e8490f", "#f18800", "#e4ce00", "#9ec417", "#13a983", "#44c1f0", "#3f60aa", "#739072", "#6F4E37")) +
    scale_color_manual(values = c("#cc340c", "#e8490f", "#f18800", "#e4ce00", "#9ec417", "#13a983", "#44c1f0", "#3f60aa", "#739072", "#6F4E37")) +
    scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
    scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) + # 修改 x 轴标签
    stat_compare_means(
      method = "t.test", paired = TRUE, comparisons = list(c("alone", "inter")),
      labels = c("alone" = "Competition-free", "inter" = "Competition")
    ) +
    labs(x = NULL, y = ylab, title = title) +
    theme_Publication() + # facet_grid(. ~ comp_year) +
    scale_linetype_manual(
      values = c("alone" = "dotted", "inter" = "solid"),
      labels = c("alone" = "Competition-free", "inter" = "Competition")
    ) +
    theme(legend.position = "none")
}


fig2a <- create_ind_trait_plot(all_indTrait_RGR_results, "Marginal_R2", "Marginal R²", "", "(a)") + 
  theme_Publication(base_size = 14) + 
  ylim(c(0, 0.28)) +
  theme(legend.position = "none")

all_indTrait_RGR_results$abs_estimate <- abs(all_indTrait_RGR_results$estimate)
fig2b <- create_ind_trait_plot(all_indTrait_RGR_results, "abs_estimate","The absolute value of the slope","","(b)") + 
  theme_Publication(base_size = 14) + 
  ylim(c(0, 0.35)) +
  theme(legend.position = "none")

combined_fig2_legend <- (fig2a | fig2b) +
  plot_layout(nrow = 1, guides = "collect") & 
  theme(legend.position = "right", 
    legend.direction = "vertical",  
    legend.box = "vertical",  
    legend.spacing.y = unit(0.5, "cm"),  
    legend.key.size = unit(0.7, "cm"), 
    axis.line = element_line(size = 1), 
    axis.ticks = element_line(size = 1),
    legend.text = element_text(size = 14),  
        legend.title = element_text(size = 14))

#combined_fig2_1

### induvidual-level
ind_tda_long <- gather(indtraits_tda, trait, values, trait.names)
ind_tda_long <- ind_tda_long[ind_tda_long$RGR > 0, ]
ind_tda_long <- na.omit(ind_tda_long)
ind_tda_long$comp_year <- factor(ind_tda_long$comp_year, levels = c("2nd", "3rd"))
ind_tda_long$Nenv <- factor(ind_tda_long$Nenv, levels = c("6", "5", "4", "9", "2", "3", "7", "8", "1"))
ind_tda_long$trait <- factor(ind_tda_long$trait, levels = trait.names)

ind_com_fig <- ggplot(ind_tda_long, aes(x = values, y = RGR, col = Compete)) +
  geom_point(alpha = 0.1, size = 1, aes(col = Compete)) +
  facet_wrap(. ~ trait,
    scales = "free", ncol = 5,
    labeller = labeller(Compete = c("alone" = "Competition-free", "inter" = "Competition"))) +
  geom_smooth(method = "lm", aes(group = Compete, fill = Compete), col = "black", alpha = 0.9, size = 0.15, se = TRUE) +
  labs(title = "(c)", x = "Trait values", y = "RGR") +
  theme_Publication() +
  scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  theme_Publication(base_size = 14) +
  theme(legend.position = "bottom",
    legend.key.size = unit(0.4, "cm"),
    legend.text = element_text(size = 15),
    legend.title = element_blank()) +
  scale_x_continuous(breaks = pretty(ind_tda_long$values, n = 5)) 

#ind_com_fig

# combined_fig2_1 <- ((fig2a / fig2b) | ind_com_fig) +plot_layout(widths = c(1, 2.5))
combined_fig2ab <- fig2a + fig2b +plot_layout(ncol = 2)
combined_fig2 <- (combined_fig2ab/ ind_com_fig) + plot_layout(heights = c(1.3, 2))
combined_fig2

topptx(combined_fig2_legend, filename = "02_output_results/figs/Fig.2_legend.pptx", height = 4, width = 8)
topptx(combined_fig2ab,filename = "02_output_results/figs/Fig.2ab.pptx", height = 4, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.png", height = 10, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.pdf", height = 10, width = 8)
ggsave(combined_fig2,filename = "02_output_results/figs/Fig.2_all_R2_slope.png", height = 10, width = 8)

################
fig_s2a <- create_ind_trait_plot(indTrait_RGR_results, "Marginal_R2", "Marginal R²",
    "", "(a)") + theme_Publication(base_size = 14) + facet_grid(. ~ comp_year) +
    theme(axis.text.x = element_blank())
indTrait_RGR_results$abs_estimate <- abs(indTrait_RGR_results$estimate)
fig_s2b <- create_ind_trait_plot(indTrait_RGR_results, "abs_estimate", "The absolute value of the slope",
    "", "(b)") + theme_Publication(base_size = 14) + facet_grid(. ~ comp_year) +
    theme(axis.text.x = element_blank())

combined_fig_s2ab <- fig_s2a + fig_s2b + plot_layout(ncol = 2, guides = "collect") &
    theme(legend.position = "bottom", legend.spacing.y = unit(0.15, "cm"), legend.key.size = unit(0.6,
        "cm"), axis.line = element_line(size = 1), axis.ticks = element_line(size = 1),
        legend.text = element_text(size = 10), legend.title = element_blank())

combined_fig_s2ab

topptx(combined_fig_s2ab, filename = "02_output_results/figs/Fig.S2_R2_slope.pptx",
    height = 5, width = 8)
ggsave(combined_fig_s2ab, filename = "02_output_results/figs/Fig.S2_R2_slope.pdf",
    height = 4.5, width = 8)

4.1 Other slopes

both_indTrait_RGR_results <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)
both_indTrait_RGR_long_re <- both_indTrait_RGR_results %>%
  pivot_longer(
    cols = c("estimate", "TLestimate", "TMestimate", "TPestimate", "TYestimate"),
    names_to = "Type",
    values_to = "estimate")%>%
  mutate(Type = case_when(
    Type == "estimate" ~ "Trait",
    Type == "TLestimate" ~ "Trait: Light",
    Type == "TMestimate" ~ "Trait: Soil moisture",
    Type == "TPestimate" ~ "Trait: Soil phosphorus",
    Type == "TYestimate" ~ "Trait: Year"))
both_indTrait_RGR_long_re$abs_estimate <- (abs(both_indTrait_RGR_long_re$estimate))

figS2b <- create_ind_trait_plot(both_indTrait_RGR_long_re[both_indTrait_RGR_long_re$comp_year=="All"&both_indTrait_RGR_long_re$Type!="Trait",], "abs_estimate","The absolute value of the slope","Treament","") + facet_wrap(.~Type,ncol=2,scale="free")+
  theme_Publication(base_size = 14) + 
  theme(legend.position = "bottom", 
    legend.spacing.y = unit(0.15, "cm"),  
    legend.key.size = unit(0.5, "cm"), 
    axis.line = element_line(size = 1), 
    axis.ticks = element_line(size = 1),
    legend.text = element_text(size = 11),  
        legend.title = element_blank())+guides(
      color = guide_legend(title = "trait"), # 显示trait的图例
      shape = "none",    # 隐藏shape的图例
      linetype = "none"  # 隐藏linetype的图例
    )

figS3 <- create_ind_trait_plot(both_indTrait_RGR_long_re[both_indTrait_RGR_long_re$Type!="Trait", ], "abs_estimate","Absolute estimates in trait—growth models","Treament","") + facet_wrap(.~Type,ncol=2,scale="free")+
  theme_Publication(base_size = 14) +facet_grid(comp_year~Type,scale="free") + 
  theme(legend.position = "bottom", 
    legend.spacing.y = unit(0.15, "cm"),  
    legend.key.size = unit(0.5, "cm"), 
    axis.line = element_line(size = 1), 
    axis.ticks = element_line(size = 1),
    legend.text = element_text(size = 11),
    axis.text.x = element_text(angle = 15),
        legend.title = element_blank())+guides(
      color = guide_legend(title = "trait"), # 显示trait的图例
      shape = "none",    # 隐藏shape的图例
      linetype = "none"  # 隐藏linetype的图例
    )+ylim(c(0,0.18))

figS3

topptx(figS3, filename = "02_output_results/figs/Fig.S3_interaction_estimate.pptx", height = 9, width = 8)
ggsave(figS3, filename = "02_output_results/figs/Fig.S3_interaction_estimate.pdf", height = 9, width = 8)

topptx(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.pptx", height = 9.5, width = 8.5)
ggsave(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.pdf", height = 9.5, width = 8.5)
ggsave(figS2b, filename = "02_output_results/figs/Fig.S2_interaction_estimate.png", height = 8.5, width = 8)

5 Figure S2 single model results

all_indTrait_RGR_results$Significance <- ifelse(all_indTrait_RGR_results$p_value <
    0.001, "***", ifelse(all_indTrait_RGR_results$p_value < 0.01, "**", ifelse(all_indTrait_RGR_results$p_value <
    0.05, "*", "")))
# Add significance stars based on p_value
indTrait_RGR_results$Significance <- ifelse(indTrait_RGR_results$p_value < 0.001,
    "***", ifelse(indTrait_RGR_results$p_value < 0.01, "**", ifelse(indTrait_RGR_results$p_value <
        0.05, "*", "")))
indTrait_RGR_results1 <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)

fig_s2a <- ggplot(indTrait_RGR_results1, aes(x = reorder(trait, Marginal_R2), y = Marginal_R2,
    fill = Compete)) + geom_bar(stat = "identity", position = "dodge") + facet_grid(. ~
    comp_year) + theme_Publication() + coord_flip() + labs(title = "(a)", x = "Trait",
    y = expression("Marginal" ~ R^2, "of trait-growth models")) + scale_color_manual(values = c("#6B8A7A",
    "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free", inter = "Competition")) +
    scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free",
        inter = "Competition")) + theme(axis.text.x = element_text(angle = 20, hjust = 1))

fig_s2b <- ggplot(indTrait_RGR_results1, aes(x = reorder(trait, Marginal_R2), y = estimate,
    fill = Compete)) + geom_bar(stat = "identity", position = "dodge") + facet_grid(. ~
    comp_year) + theme_Publication() + coord_flip() + labs(title = "(b)", x = "",
    y = "Estimate of trait in growth models") + scale_color_manual(values = c("#6B8A7A",
    "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free", inter = "Competition")) +
    scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), name = "Treatment", labels = c(alone = "Competition-free",
        inter = "Competition")) + theme(axis.text.x = element_text(angle = 20, hjust = 1))

combined_figS3 <- fig_s2a + fig_s2b + plot_layout(ncol = 2, guides = "collect") &
    theme(legend.position = "bottom", legend.key.size = unit(0.4, "cm"), legend.text = element_text(size = 12),
        legend.title = element_text(size = 12))
combined_figS3

topptx(combined_figS3, filename = "02_output_results/figs/fig_s2_r2_slopes_onlytrait_models.pptx",
    height = 5, width = 8)
ggsave(combined_figS3, filename = "02_output_results/figs/fig_s3_r2_slopes_onlytrait_models.png",
    height = 5, width = 8)

6 simple linear regression

ind_com_fig + facet_wrap(comp_year ~ trait, scales = "free", ncol = 5, labeller = labeller(Compete = c(alone = "Competition-free",
    inter = "Competition")))

7 Figure 3 trait variability and RGR

fric_indtrait <- ggplot(all_indTrait_RGR_results, aes(x = Compete, y = Fric, group = Compete)) +
  geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete), width = 0.5, alpha = 0.5) +
  scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  stat_compare_means(
    method = "t.test", paired = TRUE,
    comparisons = list(c("alone", "inter")), labels = c("alone" = "Competition-free", "inter" = "Competition")
  ) +
  labs(x = NULL, y = "Trait variability", title = "(a)") +
  theme_Publication() +
  scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  theme(legend.position = "none")#+ ylim(c(2.8, 5))

rgr_indtrait <- ggplot(indtraits_tda, aes(x = Compete, y = RGR, group = Compete)) +
  geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete), width = 0.5, alpha = 0.5) +
  scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  scale_shape_manual(values = c(1, 19), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  stat_compare_means(
    method = "t.test",
    comparisons = list(c("alone", "inter")), labels = c("alone" = "Competition-free", "inter" = "Competition")
  ) +
  labs(x = NULL, y = "Relative growth rate (RGR)", title = "(b)") +
  theme_Publication() +
  scale_x_discrete(labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  theme(legend.position = "none")#+ ylim(c(0, 1.2))

combined.fig3ab <- fric_indtrait + rgr_indtrait + plot_layout(ncol = 2, guides = "collect")&
  theme(legend.position = "bottom",
    legend.key.size = unit(0.7, "cm"),
    legend.text = element_text(size = 13),
    legend.title = element_blank())

combined.fig3ab

topptx(combined.fig3ab, filename = "02_output_results/figs/Fig.3_TV_RGR.pptx", height = 4, width = 7)
ggsave(combined.fig3ab, filename = "02_output_results/figs/Fig.3_TV_RGR.pdf", height = 4, width = 7)

itv_indtrait <- ggplot(indTrait_RGR_results, aes(x = comp_year, y = Fric, fill = Compete)) +
  geom_boxplot(position = position_dodge(width = 0.7), width = 0.5, alpha = 0.5) +  
  scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  stat_compare_means(
    method = "t.test", paired = TRUE,
    comparisons = list(c("alone", "inter")), 
    label = "p.signif") +
  labs(x = "Year under competition", y = "Trait variability", title = "(a)") +
  theme_Publication() +
  theme(axis.text.x = element_text(angle = 0, hjust = 1),  
        legend.position = "bottom", 
        legend.key.size = unit(0.4, "cm"),
        legend.text = element_text(size = 10),  
        legend.title = element_text(size = 12)) #+ ylim(c(2.8, 5))

rgr_indtrait1 <- ggplot(indtraits_tda, aes(x = comp_year, y = RGR, fill = Compete)) +
  geom_boxplot(position = position_dodge(width = 0.8), width = 0.5, alpha = 0.5) +  # 增加 position_dodge 的 width
  scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition")) +
  stat_compare_means(
    method = "t.test",
    comparisons = list(c("alone", "inter")), 
    label = "p.signif" 
  ) +
  labs(x = "Year under competition", y = "Relative growth rate (RGR)", title = "(b)") +
  theme_Publication() +
  theme(axis.text.x = element_text(angle = 0, hjust = 1)) 

combind_figS3_1 <- itv_indtrait + rgr_indtrait1 + plot_layout(ncol = 2, guides = "collect")&
  theme(legend.position = "bottom",
    legend.key.size = unit(0.7, "cm"),
    legend.text = element_text(size = 13),
    legend.title = element_blank())

combind_figS3_1

topptx(combind_figS3_1,filename = "02_output_results/figs/Fig.s3_TV_RGR.pptx", height = 4.5, width = 8)
ggsave(combind_figS3_1,filename = "02_output_results/figs/Fig.s3_TV_RGR.pdf", height = 4.5, width = 8)

8 Figure X. all diversity indexs

Fric_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = com_ITV, group = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
        width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
    labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
    19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
    paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
        inter = "Competition")) + labs(x = "Year under competition", y = "Fric",
    title = "(a)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank())  #+ ylim(c(0.05, 0.6))

Fdiv_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = Fdiv, group = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
        width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
    labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
    19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
    paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
        inter = "Competition")) + labs(x = "Year under competition", y = "FDiv",
    title = "(b)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank())  #+ ylim(c(0.05, 0.6))

indTrait_RGR_results1 <- rbind(all_indTrait_RGR_results, indTrait_RGR_results)
spitv_indtrait <- ggplot(indTrait_RGR_results1, aes(x = comp_year, y = sp_ITV, fill = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.7), width = 0.5, alpha = 0.5) +
    scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c(alone = "Competition-free",
        inter = "Competition")) + stat_compare_means(method = "t.test", paired = TRUE,
    comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
        inter = "Competition")) + labs(x = "Year under competition", y = "ITV", title = "(c)") +
    theme_Publication() + theme(axis.text.x = element_text(angle = 0, hjust = 1),
    legend.position = "bottom", legend.key.size = unit(0.4, "cm"), legend.text = element_text(size = 10),
    legend.title = element_text(size = 12))

fig_sp_itv <- ggplot(indTrait_RGR_results, aes(x = Compete, y = sp_ITV, group = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
        width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
    labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
    19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
    paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
        inter = "Competition")) + labs(x = "Year under competition", y = "Community trait variability",
    title = "(d)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank())  #+ ylim(c(0.05, 0.6))


interdis_indtrait <- ggplot(indTrait_RGR_results, aes(x = Compete, y = com_ITV, group = Compete)) +
    geom_boxplot(position = position_dodge(width = 0.5), aes(group = Compete, fill = Compete),
        width = 0.5, alpha = 0.5) + scale_fill_manual(values = c("#6B8A7A", "#DAC0A3"),
    labels = c(alone = "Competition-free", inter = "Competition")) + scale_shape_manual(values = c(1,
    19), labels = c(alone = "Competition-free", inter = "Competition")) + stat_compare_means(method = "t.test",
    paired = TRUE, comparisons = list(c("alone", "inter")), labels = c(alone = "Competition-free",
        inter = "Competition")) + labs(x = "Year under competition", y = "Community trait variability",
    title = "(d)") + theme_Publication() + facet_grid(. ~ comp_year) + theme(axis.text.x = element_blank())  #+ ylim(c(0.05, 0.6))

Fric_indtrait + Fdiv_indtrait + spitv_indtrait + interdis_indtrait + plot_layout(ncol = 2)

9 Figure 3 Mixed-effect model

# Add the 'growth_year' column based on the 'comp_year' column
indTrait_RGR_results$growth_year <- ifelse(indTrait_RGR_results$comp_year == "2nd", 2,
  ifelse(indTrait_RGR_results$comp_year == "3rd", 3, 1))

# View the updated data structure
indTrait_RGR_results$growth_year <- scale(indTrait_RGR_results$growth_year)

ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + (1 | trait),
  data = all_indTrait_RGR_results,
  family = beta_family()
)

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          
## Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##    -86.5    -80.5     49.2    -98.5       14 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.06429  0.2536  
## Number of obs: 20, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  575 
## 
## Conditional model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -2.4400     0.2016 -12.106  < 2e-16 ***
## com_ITV                1.4608     0.7543   1.937   0.0528 .  
## RGR_mean               1.0655     0.1555   6.850 7.36e-12 ***
## com_ITV:Competeinter  -0.6760     0.5728  -1.180   0.2380    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Load the glmmTMB package for advanced mixed-effects models
library(glmmTMB)
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
  Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + growth_year + (1 | trait),
  data = indTrait_RGR_results,
  family = beta_family()
)

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          
## Marginal_R2 ~ com_ITV + com_ITV:Compete + RGR_mean + growth_year +  
##     (1 | trait)
## Data: indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##   -179.0   -167.2     96.5   -193.0       33 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.1319   0.3632  
## Number of obs: 40, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  237 
## 
## Conditional model:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -2.82905    0.26345 -10.739  < 2e-16 ***
## com_ITV               1.77990    0.98745   1.803   0.0715 .  
## RGR_mean              1.13453    0.19170   5.918 3.25e-09 ***
## growth_year           0.22517    0.05346   4.212 2.53e-05 ***
## com_ITV:Competeinter -0.62508    0.68126  -0.918   0.3589    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(ind_trait_model_beta)

##            R2m       R2c
## [1,] 0.7051071 0.9282304

# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Mean growth rate", "Year", "Trait variability:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Year", "Mean growth rate", "Trait variability:Competition", "Trait variability"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
  xlab(NULL) +
  ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
  ggtitle("(c)") +
  theme_Publication(base_family = "sans") +
  coord_flip() +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 3.55, xmax = 4.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 2.55, xmax = 3.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 2.45), fill = "grey", color = NA, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)

fig_traitvar

# itv_indtrait + rgr_indtrait + fig_traitvar + plot_layout(ncol = 2, nrow = 2, design = "ABCC")
topptx(fig_traitvar,
  filename = "02_output_results/figs/Fig3.2_mixed_effect_model.pptx", height = 5, width = 7
)

library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = (com_ITV), y = Marginal_R2, color = Compete)) +
  geom_point(alpha = 1, size = 2) + # Add points with transparency for better visualization
  geom_smooth(method = "lm", size = 1, alpha = 0.5, fill="grey70") +
  geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
  labs(x = "ITV", y = "Marginal R²", color = "Compete") + # Axis and legend labels
  theme_Publication() +
  facet_grid(. ~ comp_year) +
  stat_poly_eq(
    method = "lm",
    aes(label = paste(after_stat(rr.label),
      after_stat(p.value.label),
      sep = "*\", \"*"
    )), size = 3
  ) + 
  theme(legend.position = "bottom",
    legend.key.size = unit(0.4, "cm"),
    legend.text = element_text(size = 15),
    legend.title = element_blank()) +
  scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv

---

10 Supplementary

10.1 Figure SX Fric - Mixed-effect model

# Fit a Beta regression model using glmmTMB()
all_ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait),
  data = all_indTrait_RGR_results, family = beta_family())
summary(all_ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##    -91.0    -85.0     51.5   -103.0       14 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.03922  0.198   
## Number of obs: 20, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  588 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -3.8241     0.5980  -6.395 1.61e-10 ***
## Fric                0.9575     0.3034   3.156  0.00160 ** 
## RGR_mean           -2.2323     1.2726  -1.754  0.07942 .  
## Fric:Competeinter  -0.9836     0.3564  -2.760  0.00579 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ind_trait_model_beta <- glmmTMB(
  Marginal_R2 ~ Fric + Fric:Compete + RGR_mean +(1 | trait) + (1 | comp_year),
  data = indTrait_RGR_results,  family = beta_family())

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          Marginal_R2 ~ Fric + Fric:Compete + RGR_mean + (1 | trait) +  
##     (1 | comp_year)
## Data: indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##   -175.3   -163.5     94.6   -189.3       33 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance  Std.Dev. 
##  trait     (Intercept) 1.203e-01 3.468e-01
##  comp_year (Intercept) 5.449e-11 7.382e-06
## Number of obs: 40, groups:  trait, 10; comp_year, 2
## 
## Dispersion parameter for beta family ():  205 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.00448    0.60718  -4.948 7.49e-07 ***
## Fric               0.29915    0.17226   1.737   0.0825 .  
## RGR_mean           0.10839    0.19267   0.563   0.5737    
## Fric:Competeinter -0.32091    0.05605  -5.726 1.03e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(ind_trait_model_beta)

##            R2m       R2c
## [1,] 0.6915189 0.9108413

###################################################
ind_trait_model_beta1 <- glmmTMB(Marginal_R2 ~ Fric + Fric:Compete + (1 | trait),
  data = all_indTrait_RGR_results, family = beta_family())
summary(ind_trait_model_beta1)

##  Family: beta  ( logit )
## Formula:          Marginal_R2 ~ Fric + Fric:Compete + (1 | trait)
## Data: all_indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##    -89.8    -84.8     49.9    -99.8       15 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.05413  0.2327  
## Number of obs: 20, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  567 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.29265    0.53643  -6.138 8.35e-10 ***
## Fric               0.49223    0.15401   3.196  0.00139 ** 
## Fric:Competeinter -0.35972    0.02241 -16.053  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ind_trait_model_beta2 <- glmmTMB(Marginal_R2 ~ RGR_mean + (1 | trait),
  data = all_indTrait_RGR_results,  family = beta_family())
summary(ind_trait_model_beta2)

##  Family: beta  ( logit )
## Formula:          Marginal_R2 ~ RGR_mean + (1 | trait)
## Data: all_indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##      -86      -82       47      -94       16 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.06781  0.2604  
## Number of obs: 20, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  409 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.17922    0.09136  -23.85   <2e-16 ***
## RGR_mean     1.21163    0.07922   15.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Extract coefficients and confidence intervals
coef_ind_cr1 <- broom.mixed::tidy(ind_trait_model_beta1, effects = "fixed", conf.int = TRUE)
coef_ind_cr1 <- coef_ind_cr1[coef_ind_cr1$term != "(Intercept)", ]
coef_ind_cr2 <- broom.mixed::tidy(ind_trait_model_beta2, effects = "fixed", conf.int = TRUE)
coef_ind_cr2 <- coef_ind_cr2[coef_ind_cr2$term != "(Intercept)", ]
coef_ind_cr <- rbind(coef_ind_cr1, coef_ind_cr2)
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Trait variability:Competition","Mean growth rate")

# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Mean growth rate","Trait variability:Competition", "Trait variability"))
fig3c <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
  xlab(NULL) +
  ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
  ggtitle("(c)") +
  theme_Publication(base_family = "sans") +
  coord_flip() +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 2.55, xmax = 3.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "grey60", color = NA, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)

fig3c

combined.fig3 <- (fric_indtrait|rgr_indtrait)/fig3c + 
  plot_layout(heights = c(1, 0.75))

combined.fig3

topptx(combined.fig3ab, filename = "02_output_results/figs/Fig.3ab_TV_RGR.pptx", height = 4, width = 7)
ggsave(combined.fig3ab, filename = "02_output_results/figs/Fig.3ab_TV_RGR.pdf", height = 4, width = 7)

topptx(fig3c, filename = "02_output_results/figs/Fig.3c.pptx", height = 3, width = 6)
ggsave(fig3c, filename = "02_output_results/figs/Fig.3c.pdf", height = 3, width = 6)


library(ggrepel)
figS3c <- ggplot(both_indTrait_RGR_results, aes(x = Fric, y = Marginal_R2, color = Compete)) +
   geom_point(alpha = 1, size = 2) + # Add points with transparency for better visualization
  geom_smooth(method = "lm", size = 1, alpha = 0.5, fill="grey70",se=F) +
  geom_text_repel(aes(label = trait), size = 3.6, max.overlaps = 10) +
  labs(x = "Trait variability", y = "Marginal R²", color = "Compete") + # Axis and legend labels
  facet_grid(. ~ comp_year,scale="free") +
  stat_poly_eq(
    method = "lm",
    aes(label = paste(after_stat(rr.label),
      after_stat(p.value.label),
      sep = "*\", \"*"
    )), size = 3
  ) + theme_Publication(base_size = 14) +
  theme(legend.position = "bottom",
    legend.key.size = unit(0.4, "cm"),
    legend.text = element_text(size = 13),
    legend.title = element_blank()) + ylim(c(0,0.34))+
  scale_color_manual(values = c("#6B8A7A", "#DAC0A3"), labels = c("alone" = "Competition-free", "inter" = "Competition"))

figS3c

topptx(figS3c, filename = "02_output_results/figs/Fig.S3c.pptx", height = 4, width = 8)
ggsave(figS3c, filename = "02_output_results/figs/Fig.S3c.pdf", height = 4, width = 8)

10.2 Figure SX sp_slope_cv - Mixed-effect model

ind_trait_model_beta <- glmmTMB(Marginal_R2 ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait),
  data = all_indTrait_RGR_results,
  family = beta_family()
)

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          Marginal_R2 ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait)
## Data: all_indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##    -80.8    -75.8     45.4    -90.8       15 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  trait  (Intercept) 0.06259  0.2502  
## Number of obs: 20, groups:  trait, 10
## 
## Dispersion parameter for beta family ():  294 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -2.2492     0.2641  -8.517   <2e-16 ***
## sp_slope_cv                0.8164     0.3130   2.608   0.0091 ** 
## sp_slope_cv:Competeinter  -1.5147     0.1241 -12.202   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- lmer(
  estimate ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait) + (1 | comp_year),
  data = indTrait_RGR_results
)

summary(ind_trait_model_beta)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: estimate ~ sp_slope_cv + sp_slope_cv:Compete + (1 | trait) +  
##     (1 | comp_year)
##    Data: indTrait_RGR_results
## 
## REML criterion at convergence: -70.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3867 -0.3881 -0.1127  0.2269  1.9702 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  trait     (Intercept) 0.014616 0.12090 
##  comp_year (Intercept) 0.000000 0.00000 
##  Residual              0.003786 0.06153 
## Number of obs: 40, groups:  trait, 10; comp_year, 2
## 
## Fixed effects:
##                          Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)              -0.08023    0.06320 32.39836  -1.269  0.21335   
## sp_slope_cv               0.19409    0.06682 29.35917   2.905  0.00692 **
## sp_slope_cv:Competeinter -0.06418    0.02484 28.14070  -2.583  0.01527 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) sp_sl_
## sp_slope_cv -0.767       
## sp_slp_cv:C  0.142 -0.366
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

r.squaredGLMM(ind_trait_model_beta)

##             R2m       R2c
## [1,] 0.06783866 0.8082099

# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("sp_slope_cv", "sp_slope_cv:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("sp_slope_cv:Competition", "sp_slope_cv"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
  xlab(NULL) +
  ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
  ggtitle("(c)") +
  theme_Publication(base_family = "sans") +
  coord_flip() +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)

fig_traitvar

library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = sp_slope_cv, y = estimate, color = Compete)) +
  geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
  geom_smooth(method = "lm", size = 1, alpha = 0.5) +
  geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
  labs(x = "Estimate", y = "Marginal R²", color = "Compete") + # Axis and legend labels
  theme_Publication() +
  facet_grid(. ~ comp_year) +
  stat_poly_eq(
    method = "lm",
    aes(label = paste(after_stat(rr.label),
      after_stat(p.value.label),
      sep = "*\", \"*"
    )), size = 3
  ) +
  scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv

10.3 Figure SX Fdiv - Mixed-effect model

# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
  Marginal_R2 ~ Fdiv + Fdiv:Compete + (1 | trait) + (1 | comp_year),
  data = indTrait_RGR_results,
  family = beta_family()
)

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          
## Marginal_R2 ~ Fdiv + Fdiv:Compete + (1 | trait) + (1 | comp_year)
## Data: indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##   -179.8   -169.7     95.9   -191.8       34 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance  Std.Dev. 
##  trait     (Intercept) 7.603e-02 2.757e-01
##  comp_year (Intercept) 5.473e-11 7.398e-06
## Number of obs: 40, groups:  trait, 10; comp_year, 2
## 
## Dispersion parameter for beta family ():  196 
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -3.8551     0.8240  -4.678 2.89e-06 ***
## Fdiv                3.5394     1.5105   2.343   0.0191 *  
## Fdiv:Competeinter  -2.1671     0.1722 -12.587  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(ind_trait_model_beta)

##            R2m       R2c
## [1,] 0.7587212 0.9030558

# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Fdiv", "Fdiv:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Fdiv:Competition", "Fdiv"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
  xlab(NULL) +
  ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
  ggtitle("(c)") +
  theme_Publication(base_family = "sans") +
  coord_flip() +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)

fig_traitvar

library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = Fdiv, y = Marginal_R2, color = Compete)) +
  geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
  geom_smooth(method = "lm", size = 1, alpha = 0.5) +
  geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
  labs(x = "Fdiv", y = "Marginal R²", color = "Compete") + # Axis and legend labels
  theme_Publication() +
  facet_grid(. ~ comp_year) +
  stat_poly_eq(
    method = "lm",
    aes(label = paste(after_stat(rr.label),
      after_stat(p.value.label),
      sep = "*\", \"*"
    )), size = 3
  ) +
  scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv

10.4 Figure 2 sp_ITV Mixed-effect model

# Load the glmmTMB package for advanced mixed-effects models
library(glmmTMB)
# Fit a Beta regression model using glmmTMB()
ind_trait_model_beta <- glmmTMB(
  Marginal_R2 ~ sp_ITV + sp_ITV:Compete + (1 | trait) + (1 | comp_year),
  data = indTrait_RGR_results,
  family = beta_family()
)

summary(ind_trait_model_beta)

##  Family: beta  ( logit )
## Formula:          
## Marginal_R2 ~ sp_ITV + sp_ITV:Compete + (1 | trait) + (1 | comp_year)
## Data: indTrait_RGR_results
## 
##      AIC      BIC   logLik deviance df.resid 
##   -161.3   -151.1     86.6   -173.3       34 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance  Std.Dev. 
##  trait     (Intercept) 5.576e-02 2.361e-01
##  comp_year (Intercept) 8.042e-11 8.968e-06
## Number of obs: 40, groups:  trait, 10; comp_year, 2
## 
## Dispersion parameter for beta family ():  107 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.66890    0.46022  -3.626 0.000288 ***
## sp_ITV              -0.10419    0.17661  -0.590 0.555245    
## sp_ITV:Competeinter -0.39555    0.04996  -7.918 2.42e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r.squaredGLMM(ind_trait_model_beta)

##            R2m       R2c
## [1,] 0.7593817 0.8508415

# Extract coefficients and confidence intervals
coef_ind_cr <- broom.mixed::tidy(ind_trait_model_beta, effects = "fixed", conf.int = TRUE)
coef_ind_cr <- coef_ind_cr[coef_ind_cr$term != "(Intercept)", ]
coef_ind_cr$col <- ifelse(coef_ind_cr$conf.low * coef_ind_cr$conf.high < 0, "white", "black")
colnames(coef_ind_cr)[2] <- "Parameter"
Parameter_names <- c("Trait variability", "Trait variability:Competition")
# Define parameter levels and plot
coef_ind_cr$Parameter <- factor(Parameter_names, levels = c("Year", "Mean growth rate", "Trait variability:Competition", "Trait variability"))
fig_traitvar <- ggplot(coef_ind_cr, aes(x = Parameter, y = estimate)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.4) +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3) +
  xlab(NULL) +
  ylab(expression("Effect size of " * italic("Marginal R"^2) * " in trait-growth")) +
  ggtitle("(c)") +
  theme_Publication(base_family = "sans") +
  coord_flip() +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 1.55, xmax = 2.45), fill = "#6B8A7A", color = NA, alpha = 0.1) +
  geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "#DAC0A3", color = NA, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_point(shape = 21, fill = coef_ind_cr$col, lwd = 4, size = 3)

fig_traitvar

# itv_indtrait + rgr_indtrait + fig_traitvar + plot_layout(ncol = 2, nrow = 2, design = "ABCC")
topptx(fig_traitvar,
  filename = "02_output_results/figs/Fig3.2_mixed_effect_model.pptx", height = 5, width = 7
)

library(ggrepel)
r_square_itv <- ggplot(indTrait_RGR_results, aes(x = (com_ITV), y = Marginal_R2, color = Compete)) +
  geom_point(alpha = 0.7, size = 2) + # Add points with transparency for better visualization
  geom_smooth(method = "lm", size = 1, alpha = 0.5) +
  geom_text_repel(aes(label = trait), size = 3, max.overlaps = 10) +
  labs(x = "ITV", y = "Marginal R²", color = "Compete") + # Axis and legend labels
  theme_Publication() +
  facet_grid(. ~ comp_year) +
  stat_poly_eq(
    method = "lm",
    aes(label = paste(after_stat(rr.label),
      after_stat(p.value.label),
      sep = "*\", \"*"
    )), size = 3
  ) +
  scale_color_manual(values = c("#609966", "#FFBF61"), labels = c("alone" = "Competition-free", "inter" = "Competition"))
r_square_itv

10.5 Figure SX spITV and comITV

ind_trait_model2 <- lmerTest::lmer(Marginal_R2 ~ RGR_mean + com_ITV + sp_ITV + com_ITV:Compete +
    sp_ITV:Compete + (1 | comp_year) + (1 | trait), data = indTrait_RGR_results)
tab_model(ind_trait_model2)

	Marginal_R2
Predictors	Estimates	CI	p
(Intercept)	0.09	-0.02 – 0.20	0.121
RGR mean	0.13	0.05 – 0.22	0.004
com ITV	0.15	-0.05 – 0.35	0.142
sp ITV	-0.02	-0.05 – 0.02	0.288
com ITV × Competeinter	-0.08	-0.23 – 0.07	0.277
sp ITV × Competeinter	0.02	-0.00 – 0.05	0.084
Random Effects
σ2	0.00
τ00 trait	0.00
τ00 comp_year	0.00
ICC	0.77
N comp_year	2
N trait	10
Observations	40
Marginal R2 / Conditional R2	0.501 / 0.885

# Figure
coef_ind_cr <- rbind(effectsize(ind_trait_model2)[2:6, ])
colnames(coef_ind_cr)[2:3] <- c("Std_Coefficient", "CI")
coef_ind_cr$col <- "black"
coef_ind_cr$col[coef_ind_cr$CI_low * coef_ind_cr$CI_high < 0] <- "white"
coef_ind_cr$Parameter <- factor(c("Mean RGR", "Community trait variability", "Species trait variability",
    "Community trait variability:Competition", "Species trait variability:Competition"),
    levels = c("Mean RGR", "Community trait variability:Competition", "Species trait variability:Competition",
        "Community trait variability", "Species trait variability"))

fig_traitvar_RGR <- ggplot(coef_ind_cr, aes(x = Parameter, y = Std_Coefficient)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") + geom_errorbar(aes(ymin = CI_low,
    ymax = CI_high), width = 0.4) + geom_point(shape = 21, fill = coef_ind_cr$col,
    lwd = 4, size = 3) + xlab(NULL) + ylab(expression("Effect size of " * italic("Marginal R"^2) *
    " in trait-growth")) + ggtitle("") + theme_Publication(base_family = "sans") +
    coord_flip() + geom_rect(aes(ymin = -Inf, ymax = Inf, xmin = 3.55, xmax = 5.45),
    fill = "#6B8A7A", color = NA, alpha = 0.1) + geom_rect(aes(ymin = -Inf, ymax = Inf,
    xmin = 1.55, xmax = 3.45), fill = "#DAC0A3", color = NA, alpha = 0.1) + geom_rect(aes(ymin = -Inf,
    ymax = Inf, xmin = 0.55, xmax = 1.45), fill = "grey", color = NA, alpha = 0.1) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") + geom_errorbar(aes(ymin = CI_low,
    ymax = CI_high), width = 0.4) + geom_point(shape = 21, fill = coef_ind_cr$col,
    lwd = 4, size = 3)

fig_traitvar_RGR

10.6 Figure 3-2 SEM

library(piecewiseSEM)
data.sem <- all_indTrait_RGR_results[, c("Compete", "comp_year", "com_ITV", "RGR_mean",
    "Fric", "Marginal_R2", "trait")]
data.sem$Compete <- ifelse(data.sem$Compete == "alone", 0, ifelse(data.sem$Compete ==
    "inter", 1, NA))
data.sem$Compete <- (as.numeric(data.sem$Compete))
data.sem.na <- data.sem
data.sem <- na.omit(data.sem)

sem_fit_full2 <- psem(lmer(RGR_mean ~ Compete + (1 | trait), na.action = na.omit,
    data = data.sem), lmer(com_ITV ~ Compete + (1 | trait), na.action = na.omit,
    data = data.sem), lmer(Marginal_R2 ~ Compete * com_ITV + RGR_mean + (1 | trait),
    na.action = na.omit, data = data.sem))
# summary(sem_fit_full2, standardized = TRUE, fit.measures = TRUE, rsquare =
# TRUE) Run summary (sem.sum2 <- summary(sem_fit_full2)) sem.sum2$coefficients
plot(sem_fit_full2, node_attrs = list(shape = "rectangle", color = "white", fillcolor = "white"),
    digits = 2)

########################
sem_fit_full2 <- psem(lmer(RGR_mean ~ Compete + (1 | trait), na.action = na.omit,
    data = data.sem), lmer(Fric ~ Compete + (1 | trait), na.action = na.omit, data = data.sem),
    lmer(Marginal_R2 ~ Compete * Fric + RGR_mean + (1 | trait), na.action = na.omit,
        data = data.sem))
# summary(sem_fit_full2, standardized = TRUE, fit.measures = TRUE, rsquare =
# TRUE) Run summary (sem.sum2 <- summary(sem_fit_full2)) sem.sum2$coefficients
plot(sem_fit_full2, node_attrs = list(shape = "rectangle", color = "white", fillcolor = "white"),
    digits = 2)

11 Supplementary

11.1 Figure S1

ind_sp_fig <- ggplot(ind_tda_long[ind_tda_long$trait %in% c("LTO", "LA", "SLA", "SSD",
    "RTD"), ], aes(x = (values), y = RGR, col = sp)) + geom_point(alpha = 1, size = 0.5,
    aes(col = sp)) + facet_grid(Compete ~ trait, labeller = labeller(Compete = c(alone = "Competition-free",
    inter = "Competition"))) + stat_ellipse(aes(group = sp, fill = sp), geom = "polygon",
    alpha = 0.5, color = NA, level = 0.9) + geom_smooth(method = "lm", se = F, alpha = 1,
    size = 1, aes(col = sp)) + labs(title = "", x = "Trait values", y = "Relative growth rate (RGR)") +
    stat_poly_eq(method = "lm", aes(label = paste(after_stat(rr.label), after_stat(p.value.label),
        sep = "*\", \"*")), size = 3) + theme_Publication() + scale_fill_Publication() +
    scale_colour_Publication()
ind_sp_fig